
rm(list=ls(all=TRUE))

set.seed(123)

library(systemfit)
library(plyr)
library(dplyr)
library(RColorBrewer)
library(miscTools)

setwd("replication_archive/SettlerMortality")

### Load the data
load("data_settlermortality.rda")


etavec = seq(0.98, 10.8, length.out=1000)

results <-  NULL

for(i in 1:length(etavec)){
	# create data given eta
	data$logem4_true <- NA
	data$logem4_true[data$benchmark==0] <- data$logem4[data$benchmark==0]
	data$logem4_true[data$benchmark==1] <- data$logem4[data$benchmark==1] + log(etavec[i]/4.25)

	# run the model with the simulated data: first stage
	s1 <- lm(avexpr ~ logem4_true + laborer + campaign, data = data)

	# run the model with the simulated data: IV
	m <- systemfit(logpgp95 ~ avexpr + laborer + campaign, inst = ~ logem4_true + laborer + campaign, data = data, method = "2SLS")

	# save the coefficients and standard errors of both stages
	est1 <- coef(s1)[2]
	err1 <- sqrt(vcov(s1)[2,2])

	est2 <- coef(m)[2]
	err2 <- sqrt(vcov(m)[2,2])

	results <- rbind(results, c(etavec[i], est1, err1, est2, err2))
}

colnames(results) <- c("eta", "pointest_first", "stderr_first", "pointest_second", "stderr_second")

results <- as.data.frame(results)

results$lowerci95_first <- results$pointest_first + qnorm(0.025)*results$stderr_first
results$upperci95_first <- results$pointest_first + qnorm(0.975)*results$stderr_first

results$lowerci95_second <- results$pointest_second + qnorm(0.025)*results$stderr_second
results$upperci95_second <- results$pointest_second + qnorm(0.975)*results$stderr_second




### FIGURE 6

### Panel (a): first stage
quartz(type="pdf", width=7, height=7, file="output/Bishops_stage1.pdf")
par(mar = c(4,4,0.3,0.3), mgp=c(2.5,1,0), family="CMU Serif")
plot(1:length(results$pointest_first), results$pointest_first, type="n", ylim=c(min(results$lowerci95_first), 0), xlab = expression(eta), ylab="First Stage Coefficient", xaxt="n")
polygon(c(1:length(results$pointest_first), rev(1:length(results$pointest_first))), c(results$lowerci95_first, rev(results$upperci95_first)), col="grey", border=NA)
points(1:length(results$pointest_first), results$pointest_first, type="l", lwd=3)

axis(1, at=seq(1, 1001, length.out=11), labels = round(seq(0.98, 10.8, length.out=11), 2), las=2)
abline(h=0, col = "black")
lines(c(333,333), c(results$lowerci95_first[333], results$upperci95_first[333]), lwd=3)  # 333 is position of 4.25
points(333, results$pointest_first[333], pch=16)
dev.off()



### Panel (b): second stage
quartz(type="pdf", width=7, height=7, file="output/Bishops_stage2.pdf")
par(mar = c(4,4,0.3,0.3), mgp=c(2.5,1,0), family="CMU Serif")
plot(1:length(results$pointest_second), results$pointest_second, type="n", ylim=c(0, max(results$upperci95_second)), xlab = expression(eta), ylab="Second Stage Coefficient", xaxt="n")
polygon(c(1:length(results$pointest_second), rev(1:length(results$pointest_second))), c(results$lowerci95_second, rev(results$upperci95_second)), col="grey", border=NA)
points(1:length(results$pointest_second), results$pointest_second, type="l", lwd=3)

axis(1, at=seq(1, 1001, length.out=11), labels = round(seq(0.98, 10.8, length.out=11), 2), las=2)
abline(h=0, col = "black")
lines(c(333,333), c(results$lowerci95_second[333], results$upperci95_second[333]), lwd=3)
points(333, results$pointest_second[333], pch=16)
dev.off()





### FIGURE 5

highbel = 10.8
lowbel = .98
x = seq(lowbel, highbel, .01)/4.25
z = quantile(x,seq(0, 1, .001))
beliefs = log(z)


data1 = data$logem4 + data$benchmark*beliefs[1] 
data2 = data$logem4 + data$benchmark*beliefs[251] 
data3 = data$logem4 + data$benchmark*beliefs[501]
data4 = data$logem4 + data$benchmark*beliefs[751]
data5 = data$logem4 + data$benchmark*beliefs[1001]


quartz(type="pdf", width=7, height=7, file="output/Bishops_Data1.pdf")
par(mar = c(4,4,2,2), mgp=c(2.5,1,0), family="CMU Serif")
plot(data$logem4[data$benchmark==1], data$logem4[data$benchmark==1], pch=16, xlim=c(min(data$logem4), max(data$logem4)), ylim=c(min(data5), max(data5)), col="grey", xlab="Log Mortality", ylab="Simulated Log Mortality")
points(data$logem4[data$benchmark==0], data$logem4[data$benchmark==0], pch=16)
points(data$logem4[data$benchmark==1], data1[data$benchmark==1], pch=16, col="red")
points(data$logem4[data$benchmark==1], data2[data$benchmark==1], pch=16, col="orange")
points(data$logem4[data$benchmark==1], data3[data$benchmark==1], pch=16, col="green")
points(data$logem4[data$benchmark==1], data4[data$benchmark==1], pch=16, col="darkgreen")
points(data$logem4[data$benchmark==1], data5[data$benchmark==1], pch=16, col="blue")
text(5.45, 3.65, "0.98", col="red", cex=1.5)
text(5.45, 4.9, "3.44", col="orange", cex=1.5)
text(3.9, 4.5, "5.89", col="green", cex=1.5)
text(3.9, 4.85, "8.35", col="darkgreen", cex=1.5)
text(3.9, 5.15, "10.8", col="blue", cex=1.5)
dev.off()

quartz(type="pdf", width=7, height=7, file="output/Bishops_Data2.pdf")
par(mar = c(4,4,2,2), mgp=c(2.5,1,0), family="CMU Serif")
plot(density(data$logem4), xlim=c(0,10), lwd=2, xlab="Log Mortality", main="")
points(density(data1), type="l", lwd=2, col="red")
points(density(data2), type="l", lwd=2, col="orange")
points(density(data3), type="l", lwd=2, col="green")
points(density(data4), type="l", lwd=2, col="darkgreen")
points(density(data5), type="l", lwd=2, col="blue")
text(2.8, 0.29, "0.98", col="red", cex=1.5)
text(3.5, 0.34, "3.44", col="orange", cex=1.5)
text(5.1, 0.44, "5.89", col="green", cex=1.5)
text(5.15, 0.42, "8.35", col="darkgreen", cex=1.5)
text(5.8, 0.4, "10.8", col="blue", cex=1.5)
dev.off()

